# ******** Regressions in "Substitution Between Money and Near-Money Assets". *******
# Created by Wenhao Li, 2021 CC.
# =================== Load Packages ===================
remove(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
source("Main Analysis.R")

# =================== Figure 8: KVJ Deposits, MMF Quantity ====================================
pdf( paste0(FigurePath, "KVJ_deposits_volume.pdf"), width=5, height = 4.5 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB[Year>=1970]$YearMon, TB[ Year>=1970 , list( KVJ_net_deposits, MMF_to_GDP, Deposits_GDP ) ], LegendNames = c("KVJ Net Deposits/GDP", "MMF/GDP", "Deposits/GDP"), NA_elimination = FALSE )
dev.off()

pdf( paste0(FigurePath, "KVJ_deposits_spread.pdf"), width=5, height = 4.5 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB[Year>=1974]$YearMon, TB[ Year>=1974 , list( P2_P1_spread, P2_CP_MMF, P2_CP_Tbill, CP_deposit_spread_projection ) ], LegendNames = c("P2CP--P1CP Spread", "P2CP--MMF Spread", "P2CP--T-bill Spread", "P2CP--Deposits Spread"), NA_elimination = FALSE )
dev.off()



# =================== TABLE 12: First Stage Regressions for Table 2 ====================================
# First Stage Regressions
index = is.element(1:nrow(TB), seq(3,nrow(TB),by=3)) & !is.na(TB$LiqPrem_Log) & !is.na(TB$VIX) & !is.na(TB$FFR)  # The second one picks last month in each quarter
regs_FS=list()
regs_FS[[1]] = lm(  log(DebtToDeposits)  ~  log(DebtToGDP_raw) + DebtToGDP_raw  ,  data=TB[index]  )
regs_FS[[2]] = lm(  log(Debt_to_KVJ_Deposits)  ~ log(DebtToGDP_raw) + DebtToGDP_raw,  data=TB[index]  )
stargazer( c(regs_FS ), type="text" )
stargazer( regs_FS, omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
           omit.stat =  c("adj.rsq", "ser"),  dep.var.labels.include=FALSE, 
           covariate.labels = c( "$ \\log(\\frac{\\text{Total Tsy}}{\\text{GDP}} $)", "$ \\frac{\\text{Total Tsy}}{\\text{GDP}} $" ),
           column.labels  = c("$ \\log(\\frac{\\text{Net Tsy}}{\\text{Deposits}}) $", "$ \\log(\\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}}) $")
)

# =================== TABLE 13: First Stage Regressions for Table 3 ====================================
regs_FS = list()
regs_FS[[1]] = lm( paste0("TbillToGDP_LogDiff~" ,IV_formula) , data = TB_Diff)
regs_FS[[2]] = lm( paste0("TbillToGDP_LogDiff_Lag~" ,IV_formula) , data = TB_Diff)
regs_FS[[3]] = lm( paste0("TbillToDeposits_LogDiff~" ,IV_formula) , data = TB_Diff)
regs_FS[[4]] = lm( paste0("TbillToDeposits_LogDiff_Lag~" ,IV_formula) , data = TB_Diff)
regs_FS[[5]] = lm( paste0("FFR_Diff~" ,IV_formula) , data = TB_Diff)
stargazer( regs_FS, se=RobustSE(regs_FS), omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
           omit.stat =  c("adj.rsq", "ser"),  dep.var.labels.include=FALSE,  no.space = TRUE,
           column.labels  = c("$ \\Delta \\log(\\frac{\\text{T-bill}_t}{\\text{GDP}_t}) $", "$ \\Delta \\log(\\frac{\\text{T-bill}_{t-1}}{\\text{GDP}_{t-1}}) $", "$ \\Delta \\log( \\frac{\\text{T-bill}_t}{\\text{Deposit}_t} )$",
                              "$ \\Delta \\log( \\frac{\\text{T-bill}_{t-1}}{\\text{Deposit}_{t-1}} )$", "$ \\Delta \\text{FFR}_t  $")
)
stargazer(regs_FS, type="text")

# =================== TABLE 14: Concatenated Deposits Spread ====================================
# See file "Main Analysis_concatenated_spread.R".  Run it separately to get Table 14. 

# =================== TABLE 15: Vary the projection of deposits spread ================
delta_actual = 0.339
delta_vec = c(0.3, delta_actual, 0.4, 0.5, 0.6)
index = TB_preGMM$Year>=1990
for( deposit_option in 1 ){
  regs = list()
  # deposit_option=1:  deposits total
  # deposit_option=2:  KVJ deposits
  regs = list()
  p_value_vec = c()
  R2_vec = c()
  for( iter in 1:length(delta_vec) ){
    delta = delta_vec[iter]
    TB_delta = TB_preGMM
    TB_delta[ index, Deposit_Spread := Deposits_Spread_projection ]
    TB_delta[!index, Deposit_Spread := delta*FFR ]
    if( deposit_option==1 ){
      TBs = list( TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ])
    }else{
      TBs = list( TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ])
    }
    result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE)
    regs[[iter]] = result_with_IV$regs[[1]]
    p_value_vec = c(p_value_vec,result_with_IV$p_value_vec)
    R2_vec = c( R2_vec, result_with_IV$R2_vec )
  }
  stargazer( regs , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="15pt",
             covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines =
               list( c("p-value of J-test",  p_value_vec ),
                     c( "Variations explained", paste0(round( R2_vec*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
             column.labels  = paste0("$\\delta=$", delta_vec) ,
             dep.var.caption = "Measurement of B/D" )
}

# show the plot of R2 and p-value against the delta value
delta_vec = c(0.1, 0.2, 0.3, delta_actual, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
p_value_vec = rep(0, length(delta_vec))
R2_vec = rep(0, length(delta_vec))

for( iter in 1:length(delta_vec) ){
  delta = delta_vec[iter]
  TB_delta = TB_preGMM
  TB_delta[ index, Deposit_Spread := Deposits_Spread_projection ]
  TB_delta[!index, Deposit_Spread := delta*FFR ]
  TBs = list( TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ])
  result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE)
  regs[[iter]] = result_with_IV$regs[[1]]
  p_value_vec[iter] = c(result_with_IV$p_value_vec)
  R2_vec[iter] = c(result_with_IV$R2_vec )
}


# =================== Figure 9: Model Fit under Different Deposit Spread Projection Coefficients ================

pdf( paste0(FigurePath,"delta_and_GMM_IV_results.pdf"), width = 8, height = 4 )
par(mar = c(3.5, 3.5, 1, 1), mfrow=c(1,2))
MyPlot(delta_vec,p_value_vec, xlab=expression(delta), ylab=expression(p-value))
abline( v=delta_actual )
MyPlot(delta_vec,R2_vec, xlab=expression(delta), ylab=expression(R^2))
abline( v=delta_actual )
dev.off()


# ===================== TABLE: Other Measures of the Liquidity Premium ==========================
# What about other measures? 
TBs = list()
TBs = list( TB_preGMM[   , list( CD_Tsy_3M, Deposit_Spread, DebtToGDP, DepositsCore_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[   , list( NoteTBillSpread, Deposit_Spread, DebtToGDP, DepositsCore_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[   , list( GSW_Tsy_3M, Deposit_Spread, DebtToGDP, DepositsCore_GDP, VIX , DebtToGDP_raw) ]
)
result = core_GMM_fun(TBs, use_debt_GDP_IV=TRUE, beta_init=c(0.6,0.01), type="twoStep" )
stargazer( result$regs , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines = 
             list( c("p-value of J-test",  result$p_value_vec ), 
                   c( "Variations explained", paste0(round(result$R2_vec*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = c("CD-Tbill 3M", "Note-Bill Spread", "GSW-Tbill 3M") ,
           dep.var.caption = "Treasury Liquidity Premium Measures" )



# ===================== TABLE 17: GMM Estimation with Alternative Measures  ==========================
g_single_other_lambda <- function(beta,X_matrix,prediction=FALSE)
{
  X_matrix = as.matrix(X_matrix)
  rho = beta[1]
  beta_lambda = beta[2]
  rho_lambda = beta[3]
  # The structure of X_matrix = [  liquidity_premium, Deposit_spread, bond/GDP, deposit/GDP, consumption, GDP, VIX   ] 
  LiqPrem = X_matrix[,1]
  deposits_spread = X_matrix[,2]
  bonds_GDP = X_matrix[,3]
  deposits_GDP = X_matrix[,4]
  proxy = X_matrix[,5]
  total_debt_to_GDP = X_matrix[,6]
  FFR = X_matrix[,7]
  model_generated = deposits_spread*beta_lambda*proxy^rho_lambda*( bonds_GDP/deposits_GDP )^(rho-1) 
  residual = LiqPrem - model_generated
  iv1 = proxy
  iv2 = FFR  # Important: use the same specification
  iv3 = total_debt_to_GDP
  iv4 = total_debt_to_GDP^2
  mm = cbind( residual, residual*iv1, residual*iv2, residual*iv3, residual*iv4 )
  if(!prediction){
    return(mm)
  }else{
    return( data.frame( LiqPrem=LiqPrem, liq_prem_predicted=model_generated ) )
  }
}
GMM_fun_other_lambda <- function(TBs){
  N = length(TBs)
  R2_vec = rep(0, N)
  p_value_vec = rep(0, N)
  regs = list()
  for(iter in 1:N){
    DT = remove_missing(as.data.frame(TBs[[iter]]))
    reg = gmm_customized( function(beta,X) {g_single_other_lambda(beta, X)}, DT, c(0.4,0.01,1), optfct = "nlminb" )
    result = g_single_other_lambda(reg$coefficients, DT, prediction = TRUE)
    R2_vec[iter] = R2_fun( result$liq_prem_predicted, result$LiqPrem )
    p_value_vec[iter] = round( summary(reg)$stest$test[2], digits=3  )
    regs[[iter]] = reg
  }
  return( list( regs=regs, R2_vec=R2_vec, p_value_vec=p_value_vec )  )
}
TB_preGMM = TB[  !is.na(LiqPrem) & !is.na(DebtToGDP) & !is.na(Deposits_GDP) & !is.na(VIX) & !is.na(FFR) ,  ]
TB_preGMM = TB_preGMM[Is_Quarter_End==TRUE]
TBs = list( TB_preGMM[   , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw, FFR) ],
            TB_preGMM[   , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw, FFR) ],
            TB_preGMM[   , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, credit_spread , DebtToGDP_raw, FFR) ],
            TB_preGMM[   , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, credit_spread , DebtToGDP_raw, FFR) ],
            TB_preGMM[   , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, Lvg , DebtToGDP_raw, FFR) ],
            TB_preGMM[   , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, Lvg , DebtToGDP_raw, FFR) ]
)
result = GMM_fun_other_lambda(TBs)
stargazer( result$regs , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$", "$\\kappa$" ), add.lines = 
             list( c("p-value of J-test",  result$p_value_vec ), 
                   c( "Variations explained", paste0(round( result$R2_vec*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = rep( c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),3) ,
           dep.var.caption = "Measurement of B/D" )


# =================== TABLE 18: Use T-bill quantity instead of total quantity ====================================
# Note: results are very pool using Tbill quantity.
Tbill_GMM_fun <- function(TBs, use_debt_GDP_IV=FALSE, shut_off_VIX=FALSE, beta_init=c(0.6,0.005), lower=rep(0,2), upper=rep(1,2)){
  N = length(TBs)
  R2_vec = rep(0, N)
  p_value_vec = rep(0, N)
  regs = list()
  for(iter in 1:N){
    DT = remove_missing(as.data.frame(TBs[[iter]]))
    reg = gmm_customized( function(beta,X) {g_single(beta, X, use_debt_GDP_IV = use_debt_GDP_IV, shut_off_VIX=shut_off_VIX)} , DT, beta_init, optfct = "nlminb", lower=lower, upper=upper )
    result = g_single(reg$coefficients, DT, prediction = TRUE, use_debt_GDP_IV = use_debt_GDP_IV, shut_off_VIX=shut_off_VIX)
    R2_vec[iter] = R2_fun( result$liq_prem_predicted, result$LiqPrem )
    p_value_vec[iter] = round( as.numeric(summary(reg)$stest$test[2])  , digits=3  )
    regs[[iter]] = reg
  }
  return( list( regs=regs, R2_vec=R2_vec, p_value_vec=p_value_vec )  )
}
TBs = list( TB_preGMM[   , list( LiqPrem, Deposit_Spread,  TbillToGDP , DepositsCore_GDP, VIX , DebtToGDP_raw) ], 
            TB_preGMM[  , list( LiqPrem, Deposit_Spread, TbillToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ]
)
result = Tbill_GMM_fun(TBs)
result_with_IV = Tbill_GMM_fun(TBs, use_debt_GDP_IV = TRUE)
stargazer( c(result$regs, result_with_IV$regs), omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines = 
             list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ), 
                   c( "Total Treasury IV?", "No", "No", "Yes", "Yes" ),
                   c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = rep(c( "$ \\frac{\\text{T-bills}}{\\text{Deposits}} $", "$ \\frac{\\text{T-bills}}{\\text{KVJ Deposits}} $"),2) , 
           dep.var.caption = "Measurement of B/D" )




# =================== TABLE 19: Different Estimation Methods ====================
R2_GMM_fun <- function(reg, X_matrix, g_fun){
  beta = reg$coefficients
  result = g_fun( beta, X_matrix )
  R2 = R2_fun( result$liq_prem_predicted,  result$LiqPrem )
  return(R2)
}

beta_init = c(0.4, 0.1)
rho_estimated = 0.66
beta_lambda_estimated = 0.01
g_single_comparison <- function(beta,X_matrix, prediction=FALSE, combination=1 )
{
  X_matrix = as.matrix(X_matrix)
  rho = beta[1]
  beta_lambda = beta[2]
  # The structure of X_matrix = [  liquidity_premium, Deposit_spread, bond/GDP, deposit/GDP, consumption, GDP, VIX   ] 
  LiqPrem = X_matrix[,1]
  deposits_spread = X_matrix[,2]
  bonds_GDP = X_matrix[,3]
  deposits_GDP = X_matrix[,4]
  VIX = X_matrix[,5]
  DebtToGDP_total = X_matrix[,6]
  if(!prediction){
    residual = LiqPrem - deposits_spread*beta_lambda*VIX*( bonds_GDP/deposits_GDP )^(rho-1) 
    iv1 = deposits_spread
    iv2 = VIX
    iv3 = DebtToGDP_total
    iv3_extra = DebtToGDP_total^2
    mm = cbind( residual, residual*iv1, residual*iv2, residual*iv3, residual*iv3_extra )
    if( combination==1 ){
      return( mm[,c(1,2)] )
    }else if( combination==2 ){
      return( mm[,c(1,3)] )
    }else if( combination==3 ){
      return( mm[,c(1,4,5)] )
    }else if( combination==4 ){
      return( mm[,c(1,2,3)] )
    }else if( combination==5 ){
      return( mm[,c(1,2,4,5)] )
    }else if( combination==6 ){
      return( mm[,c(1,3,4,5)] )
    }else{
      return(mm)
    }
  }else{
    lambda = beta_lambda * VIX
    liq_prem_predicted = deposits_spread * (lambda*( bonds_GDP/deposits_GDP )^(rho-1) )
    return( data.table(liq_prem_predicted, LiqPrem) )
  }
}

DT = TB_preGMM[ , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX, DebtToGDP_raw ) ]
regs_whole_sample = list()
for(i in 1:7){
  regs_whole_sample[[i]] = gmm_customized( function(beta,X) g_single_comparison(beta,X,combination=i), DT , beta_init, optfct = "nlminb", lower=c(0, 0), upper=c(1,1) )
}
results = g_single( regs_whole_sample[[7]]$coefficients, DT, prediction=TRUE )
reg_NLS = nls(LiqPrem ~ beta_lambda*VIX*(DebtToGDP/Deposits_GDP)^(rho-1)*Deposit_Spread, data = DT, start = list(rho=rho_estimated, beta_lambda = beta_lambda_estimated) , algorithm = "port" )
regs_whole_sample[[8]] = regs_whole_sample[[7]]
regs_whole_sample[[8]]$coefficients = coef(reg_NLS)
regs_whole_sample[[8]]$vcov = vcov(reg_NLS)
R2_vec = rep(0, 8)
p_values = rep(0,8)
for( i in 1:length(R2_vec) ){
  names(regs_whole_sample[[i]]$coefficients) <- c( "rho", "beta_lambda" )
  R2_vec[i] = R2_GMM_fun( regs_whole_sample[[i]], DT, function(beta,X) g_single_comparison(beta,X, prediction=TRUE, combination=i) )
  p_values[i] = summary(regs_whole_sample[[i]])$stest$test[,2]
}
p_values[1:3] = NaN
p_values = as.numeric(p_values)
stargazer(regs_whole_sample, se=RobustSE(regs_whole_sample) , 
          df = FALSE, digits = 2, column.sep.width="0pt", covariate.labels = c("$\\rho$", "$\\beta_\\lambda$"), 
          add.lines = list( c("Variation explained", round(R2_vec,digits = 2) ),  
                            c("p-value", rep("NA",3), round(p_values[4:7],digits=2), "NA" ) ),
          star.cutoffs = c(0,0,0))

stargazer(regs_whole_sample, type="text", se=RobustSE(regs_whole_sample) , covariate.labels = c("$\\rho$", "$\\beta_\\lambda$"), 
          add.lines = list( c("Variation explained", round(R2_vec,digits = 2) ),  
                            c("p-value", rep("NA",3), round(p_values[4:7],digits=2), "NA" ) ),
          star.cutoffs = c(0,0,0))



# =================== TABLE 20-21: Baseline GMM with More Lags ====================

# =================== TABLE 20: GMM Estimation with 12 Lags for HAC Standard Errors ====================

# Robustness checks with more lags
TBs = list( TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ], 
            TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ]
)
# ==== Appendix Table ====
result = core_GMM_fun(TBs, lags=12)
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE, lags=12)
stargazer( c(result$regs, result_with_IV$regs), omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines = 
             list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ), 
                   c( "Total Treasury IV?", "No", "No", "Yes", "Yes" ),
                   c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) , 
           dep.var.caption = "Measurement of B/D" )

# =================== TABLE 21: GMM Estimation of 40 Lags for HAC Standard Errors ====================
result = core_GMM_fun(TBs, lags=40)
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE, lags=40)
stargazer( c(result$regs, result_with_IV$regs), omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines = 
             list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ), 
                   c( "Total Treasury IV?", "No", "No", "Yes", "Yes" ),
                   c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) , 
           dep.var.caption = "Measurement of B/D" )



# ===================== Table 22: Money Demand Regressions with Constant lambda and mu =======================
rho = reg_nested_baseline$coefficients[1]
beta_lambda = reg_nested_baseline$coefficients[2]
lambda =  mean_fun(beta_lambda*TB$VIX / (1+beta_lambda*TB$VIX))
eta = 1
beta_mu = reg_nested_baseline$coefficients[4]
mu = mean_fun(beta_mu*TB$VIX / (1+beta_mu*TB$VIX))
# Plot the slope coefficient on money demand. 
TB[ , logsp:=log(Deposit_Spread/(1+FFR/100)) ]
TB[ , loglsp:=log(LiqPrem/(1+FFR/100)) ]
TB[ , logY:=log(GDP_Real) ] # Log real GDP
TB[ , logcu := log(CURRCIR/GDP_nominal*GDP_Real) ] # currency in circulation
TB[ , Q1 := ((1-lambda)*(Deposits_GDP*GDP_Real)^rho + lambda*(DebtToGDP*GDP_Real)^rho)^(1/rho) ]
TB[ , NB :=  KVJ_fin_sector_debt_to_GDP - Deposits_GDP]
TB[ , B_prime := ( (1-mu)*(NB*GDP_Real)^eta + mu*(DebtToGDP*GDP_Real)^eta )^(1/eta)]
TB[ , Q2 := ( (1-lambda)*(Deposits_GDP*GDP_Real)^rho + lambda*(B_prime)^rho )^(1/rho)]
TB[ , logQ :=  log(Q1)-logY ] 
TB[ , logQ2 := log(Q2)-logY ] 
TB[ , logB := log(DebtToGDP) ]   
TB[ , logm := log(CURRCIR/GDP_nominal*GDP_Real + Checking_GDP*GDP_Real)-logY ]  # Including both currency and checking in the money demand
TB[ , logMK := log(KVJ_fin_sector_debt_to_GDP*GDP_Real)-logY ]
TB[ , logMD := log(Deposits_GDP*GDP_Real)-logY ]
TB[ , logMZ := log(MZM_to_GDP*GDP_Real)-logY ]
TB[ , logvix := log(VIX) ]
TB[ , adjsp := logsp - (rho-1)*log(KVJ_fin_sector_debt_to_GDP*GDP_Real) ]
TB[ , spread := LiqPrem/(1+FFR/100) ]
TB[ , QoverY := Q1/GDP_Real ]
TB[ , Q2overY := Q2/GDP_Real ]
TB_for_reg = TB[seq(1,nrow(TB),by=3)]
regs = list()
regs[[1]] = lm( logm ~ logsp + logY , data = TB_for_reg[Year<1980 ]  )
regs[[2]] = lm( logm ~ logsp + logY , data = TB_for_reg[Year>=1980]  )
regs[[3]] = lm(logQ ~ logsp + logY + logvix, data= TB_for_reg[Year<1980] )
regs[[4]] = lm(logQ ~ logsp + logY + logvix, data= TB_for_reg[Year>=1980] )
regs[[5]] = lm(logQ2 ~ logsp + logY + logvix, data= TB_for_reg[Year<1980] )
regs[[6]] = lm(logQ2 ~ logsp + logY + logvix, data= TB_for_reg[Year>=1980] )
stargazer(regs, se=RobustSE(regs), type="text", omit.stat =  c("F", "adj.rsq", "ser"), omit="Constant" )
stargazer(regs, se=RobustSE(regs), covariate.labels = c( "$\\log( \\frac{\\delta i_t}{1+i_t})$", "$\\log(\\text{real GDP}_t)$", "$\\log(\\text{VIX}_t)$" ), 
          star.cutoffs = c(0,0,0), dep.var.labels.include = FALSE, omit.stat =  c("F", "adj.rsq", "ser"), 
          column.labels = c(rep( c("pre 1980", "post 1980"), 3 ))  )





# ===================== Table 23: Add Currency Into Money Demand Regressions =======================
regs = list()
rho = reg_baseline$coefficients[1]
beta_lambda = reg_baseline$coefficients[2]
lambda = beta_lambda*TB$VIX / (1+beta_lambda*TB$VIX)
beta_mu = reg_nested_baseline$coefficients[4]
eta = 1
mu = beta_mu*TB$VIX / (1+beta_mu*TB$VIX)
TB[ , logsp:=log(Deposit_Spread/(1+FFR/100)) ]
TB[ , loglsp:=log(LiqPrem/(1+FFR/100)) ]
TB[ , logY:=log(GDP_Real) ] # Log real GDP
TB[ , logcu := log(CURRCIR/GDP_nominal*GDP_Real) ] # currency in circulation
TB[ , Q1 := ((1-lambda)*(Deposits_GDP*GDP_Real)^rho + lambda*(DebtToGDP*GDP_Real)^rho)^(1/rho) ]
TB[ , NB :=  KVJ_fin_sector_debt_to_GDP - Deposits_GDP]
TB[ , B_prime := ( (1-mu)*(NB*GDP_Real)^eta + mu*(DebtToGDP*GDP_Real)^eta )^(1/eta)]
TB[ , Q1bar := ((1-lambda)*((Deposits_GDP+CURRCIR/GDP_nominal)*GDP_Real)^rho + lambda*(DebtToGDP*GDP_Real)^rho)^(1/rho)  ]
TB[ , Q2bar := ( (1-lambda)*((Deposits_GDP+CURRCIR/GDP_nominal)*GDP_Real)^rho + lambda*(B_prime)^rho )^(1/rho)]
TB[ , logQbar:= log(Q1bar)-logY ]
TB[ , logQ2bar:= log(Q2bar)-logY ]
TB[ , logm := log(CURRCIR/GDP_nominal*GDP_Real + Checking_GDP*GDP_Real)-logY ]  # Including both currency and checking in the money demand
TB[ , logvix := log(VIX) ]
TB_for_reg = TB[seq(1,nrow(TB),by=3)]
regs[[1]] = lm( logm ~ logsp + logY , data = TB_for_reg[Year<1980 ]  )
regs[[2]] = lm( logm ~ logsp + logY , data = TB_for_reg[Year>=1980]  )
regs[[3]] = lm(logQbar ~ logsp + logY + logvix, data= TB_for_reg[Year<1980] )
regs[[4]] = lm(logQbar ~ logsp + logY + logvix, data= TB_for_reg[Year>=1980] )
regs[[5]] = lm(logQ2bar ~ logsp + logY + logvix, data= TB_for_reg[Year<1980] )
regs[[6]] = lm(logQ2bar ~ logsp + logY + logvix, data= TB_for_reg[Year>=1980] )
stargazer(regs, se=RobustSE(regs), type="text", omit.stat =  c("F", "adj.rsq", "ser"), omit="Constant" )
stargazer(regs, covariate.labels = c( "$\\log( \\frac{\\text{deposit spread}_t}{1+i_t})$", "$\\log(\\text{real GDP}_t)$", "$\\log(\\text{VIX}_t)$" ), 
          star.cutoffs = c(0,0,0), dep.var.labels.include = FALSE, omit.stat =  c("F", "adj.rsq", "ser"), 
          column.labels = c(rep( c("pre 1980", "post 1980"), 3 ))  )

